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Abstract 

We demonstrate that variations in molecular chemisorption energy on different metals, different 
surface terminations, and different strain conditions can be accounted for by orbital-specific changes 
in the substrate electronic structure. Our density functional theory data set, spanning three metals, 
two surface terminations, and five strain states, is fit to a single model based on tight binding. 
A crucial aspect of the model is decomposition of the ci-band into contributions from the five d 
atomic orbitals. This provides a representation of the energy levels of the substrate that are directly 
relevant to the chemisorption bond, leading to accurate prediction of chemisorption trends. 
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Currently, great attention is focused on elucidating how surface modification affects sur- 
face reactivity. Recent research shows that small changes in surface electronic 
structure, induced by alloying or strain, can significantly change surface- catalyzed reaction 
rates, [y] A quantitative understanding of how changes in surface geometry and electronic 
structure affect surface reactivity will enable the design of more specific and more effec- 
tive catalysts. It has been shown that 64, the center of the transition metal (TM) <i-band 
density of states projected on the surface atoms (PDOSV, is generally predictive of trends 
in chemisorption energies (E chcm ) on TM surfaces. |?], 0, S] However, quantitative accuracy 
(model predictions accurate within 0.1 eV) is still elusive, and for several cases there is poor 
or no correlation between and E'chem- llil llll uM ^ n this paper, using CO chemisorption 
as an example, we show that more rigorous examination of the surface electronic structure 
coupled with a simple modification of current chemisorption modeling enables us to achieve 
this goal. 

We have compiled a database of DFT molecular top site (.E'chem ) and dissociative bridge 
site (-Edissoc) chemisorption energies and electronic structure measurements for CO on Pt, 
Pd, and Rh (111) and (100) surfaces. E c ^ em and E'dissoc are determined for each surface 
at the preferred theoretical lattice constants as well as under in-plane strains of ±1% and 
±2%, a range easily achievable through epitaxial mismatch. [13| Studying the response of 
chemisorption to strain as well as to different metals and facets deepens the study. Strain 
induces relatively subtle changes in E c ^ era and -Edissoc (compared to changes in metal or 
facet), so accurately accounting for large and small changes is a stringent test of a proposed 
theoretical model. Since lateral stress changes inter-planar separations, straining the systems 
also probes the interplay between in-plane and inter-plane perturbations to the surface 
geometry. 

For each metal, surface, and strain state, two values for E chem are determined, E^ em and 
£ehem- The former is the energy gain when the same chemisorption geometry is fixed over 
the relaxed bare surface for all metals and surfaces. In the latter, full ionic relaxation is 
allowed in the top two metal surface layers and all C and O ionic degrees of freedom. For 
the dissociative systems, we only determine Ef* soc , due to the known instability of C and O 
atomic adsorption at bridge sites 

We focus attention on top-site -Echem and bridge-site E? d is SO c f° r clarity. The symmetries 
at these sites provide for zero overlap between some d orbitals and the adsorbate orbitals, 
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making a clear distinction between orbital-specific and orbital-averaged models. Analysis 
of the top site also facilitates separation of the molecular chemisorption into a and 7r con- 
tributions. However, even at low-symmetry sites, the contributions of the five d orbitals to 
chemisorption are unequal, making our orbital-specific treatment more accurate in general. 

DFT calculations are performed with a generalized-gradient approximation (GGA) 
exchange-correlation functional and norm-conserving optimized pseudopotentials H| 
with the designed nonlocal method for metals. [3, [3 Pseudopotentials were designed us- 
ing the OPIUM pseudopotential package. 1 1911 All E chem values have been corrected using 
our first-principles extrapolation procedure. |2£j] Metal surfaces are modeled as slabs of five 
layers separated by vacuum, with the c(4 x 2) surface cell for (111) surfaces and the p(2 x 2) 
surface cell for (100). CO top site and C and O bridge site chemisorption are modeled at 
coverage of 0=1/4. Calculations are done, and values of -E c hem tested to be converged within 
0.02 eV, using an 8 x 8 x 1 grid of Monkhorst-Pack fc-points, reduced by symmetry where 
possible. I21I ] 

The PDOS for each orbital is constructed by projecting each atomic valence pseudo- 
wavefunction (radial wavefunction multiplied by real combination of spherical harmonics) 
of the surface atoms onto all the Kohn-Sham orbitals. Values of are then calculated as 
the first moment of each PDOS. 

To reduce PDOS contributions from neighboring surface atoms, projection is performed 
within a sphere of radius r cut centered about the surface atom of interest. Standard practice 
is to use a constant value for r cut when comparing the PDOS and associated values of 
different surfaces, and r cut = 2 a.u. is the default in some widely used DFT packages. [22, 
I23I I24I ] However, this approach leaves significant contributions from the orbitals of other 
atoms, making the calculated value of dependent on r cut , which is undesirable. 

To eliminate contamination from the orbitals of neighboring atoms, we evaluated at 
various r cut values. (When r cut < 0.5 a.u., the number of FFT grid points is too small 
to allow spherical sampling, so data are presented for r cut > 0.6 a.u.) Figure Q shows the 
variation of with r cut for Pt(lll) and Pt(100). It is apparent that the asymptotic behavior 
of €d is not reached until r cut <C 2 a.u. Furthermore, values for different surfaces, strain 
states and orbital angular momenta depend differently on r cut . To obtain accurate values, 
we fit a purely quadratic function to the data and extrapolate to r cut = 0. This procedure 
greatly reduces the contribution from the orbitals of neighboring atoms, making comparison 
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of ed for various systems more meaningful. 

While our data confirm that -E c hem qualitatively tracks with e<j for the (111) surfaces, they 
reveal shortcomings of using for modeling bonding on different facets. From Figure |2J 
we see that -E c hem values on different facets differ by 0.06-0.25 eV (0.15-0.85 eV for ^dissoc) 
even though they have the same and the same metal. Figure 121 and Table show this is 
still true for with r cut — > a.u. Therefore, neither E^ em nor E^ ssoc can be fit as a single 
function of e^, when both (111) and (100) facets are considered. 

By contrast, we find that facet dependence of the chemisorption energies can be fit as a 
single linear function of e xzyz , the band center of the d xz and d yz orbitals. Figure El shows 
that single linear regressions of data for both facets are accurate to within 0.05 eV in all 
cases. This result demonstrates that focusing on the metal orbitals involved in bonding 
simplifies the observed chemisorption behavior and enables robust modeling. 

Similarly, Table |l] shows that the response of -E cll em and -Edissoc to strain is reflected in 
^-xzyz but not 6d- The tunability of E'chem and -Edissoc through strain is two to ten times 
greater on the (111) surfaces. However, ded/ds is identical within computational precision 
for the two facets of each metal, so ded/ds is uncorrelated with dE c ^ em /ds and dEdi SS oc/ds. 
The chemisorption tunability trend is strongly reflected in de xzyz /ds, suggesting that the 
response of the d xz and d yz orbitals plays a key role in modeling top-site molecular and 
bridge-site dissociative adsorption of CO. 

TABLE I: E^ em , E^ em , E%? ssoc , e d , and e xzyz in eV for the (111) and (100) surfaces of Pt, Rh, and 
Pd. The slope of each quantity with respect to in-plane strain, (^), is also reported, in units of 
meV/%strain. (e<j and e xzyz use r cut — > a.u.) 





chem 


dE 
ds 


chem 


dE 
ds 


dissoc 


dE 
ds 


, de 
td ds 


t xzyz ds 


Pt(lll) 


-1.49 


-36 


-1.58 


-46 


1.99 - 


158 


-1.77 20 


-1.74 18 


Pt(100) 


-1.69 


-3 


-1.81 


-5 


1.03 


-65 


-1.76 20 


-1.63 7 


Rh(lll) 


-1.56 


-27 


-1.68 


-24 


0.48 - 


120 


-1.38 31 


-1.35 25 


Rh(100) 


-1.65 


-18 


-1.74 


-9 


0.11 


-74 


-1.36 29 


-1.28 14 


Pd(lll) 


-1.23 


-25 


-1.26 


-29 


1.66 


-70 


-1.58 32 


-1.54 26 


Pd(100) 


-1.32 


-5 


-1.36 


-6 


1.50 


-29 


-1.56 19 


-1.50 9 



To explain the observed trends in E c h em as the metal identity, facet, and strain state 
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are varied, we reconsider the model for CO chemisorption. A second-order perturbative 
picture of chemisorption involving interacting molecular and metal orbitals is simple and 



intuitive, and has been known for some time. [25J, |26J] Hammer, Morikawa, and N0rskov 
(HMN) achieved significant success modeling molecular interactions with solid surfaces with 
a single perturbative term involving e^, the (i-band center. |7| Guided by the trends presented 
above, we cast an orbital-specific analysis (OS) in the HMN model form, modeling top- 
site chemisorption as a perturbative interaction between molecular orbitals and the (i-band 
PDOS of each spatial orbital: 



E chem = E sp - 4 { — ~ h fS n V n 



;i - m 



+ (l + f)3 ff VA (1) 
e* 1 - e 5(T J 

where / is the idealized filling of the metal d bands, V and S are perturbation matrix 
elements and overlap integrals, respectively, labeled by symmetry, e 27r » and e 5(J are the CO 
molecular orbital eaerg.es, W J,he band center of the d„ and i v orb.tals and e„ is the 
band center of the d z 2 orbitals. |30J] As in the original HMN model, a and j3 are introduced 
as fitting parameters common to all metals, and V£ ~ j3V^ d and S n ~ —aV n . From our 
analysis of DFT orbitals, we find that S a /S n to be sensitive in the limit of desired accuracy 
to both metal identity and adsorption geometry. Our overlap analysis gives Sa/S^ is 1.182, 
1.156, and 1.200 for Pt, Rh, and Pd, respectively. E sp is found to be -0.15 eV from DFT 
calculations on Al surfaces |27( and assumed to be independent of metal identity, facet, or 
strain. For the (111) and (100) surfaces, we find that / for each of the decomposed <i-bands 
is well approximated by the idealized filling of the metal (i-bands, / = {y — 1)/10, where v 
is the valence of the metal atom. 

The corresponding conventional orbital-averaged (OA) model form is given by: 

(-2-k* — e d 



+ (1 + f)S a V a (2) 



We fit the data for both fix and rlx chemisorption systems to the Equation Q and Equa- 
tion El [31] Figure H3 shows the correlation between DFT and model values for the more 
realistic rlx chemisorption systems. We calculate the root-mean-square error (RMSE) as 



an evaluation of the model, considering the (111) and (100) data separately and combined. 
When EquationEJis used to fit the data, the RMSE is 0.051 and 0.100 eV for the (111) and 
(100) surfaces respectively, and 0.079 eV overall. When Equation [T] is used to fit the data, 
the RMSE is 0.052 and 0.051 eV for the (111) and (100) surfaces respectively, and 0.052 eV 
overall. This shows that the more sophisticated model form of Equation Q is required to 
achieve the same level of accuracy for these two surfaces. 

We now address why different levels of model sophistication are needed to achieve the 
same accuracy in predicted -E^em on the (111) and (100) surfaces. First, we consider the 
salient electronic structure differences between the (111) and (100) surface facets. The dd 
metal bonding can be decomposed by symmetry into a, ir, and 5 contributions. 28] The 
square lattice of the (100) surface allows for strong dda overlap between neighboring d x 2_ y 2 
orbitals. Our DFT data show that e x 2_ y 2 is significantly lower on (100) surfaces than on 
(111), as shown for Pt in Figure ^ Since (100) surface atoms have eight nearest neighbors 
while the (111) atoms have nine, the other d orbitals are less stable on (100) surfaces to 
varying extents, with e xz and e yz significantly higher on (100) than (111). Averaging all 
the d orbitals causes the rise in bonding-relevant e xz and e yz to be masked by the drop in 
bonding-irrelevant e x 2_ y 2, so that even though e xzyz closely tracks the increase in -E^em on 
the (100) surface relative to the (111) surface, the averaged does not. 

When tensile strain is applied to a TM surface, the weakened in-plane bonding destabilizes 
the d orbitals. shifts upward, leading to stronger -Echem- IE] This basic prediction was 
confirmed in the DFT study of Mavrikakis et al. |9] However, different d orbitals shift by 
different amounts, based on their orientations relative to the surface. 

In addition, the inter-planar spacing between the top two metal layers (r 12 ) responds 
to the strain, and this further affects substrate electronic structure, again in an orbital- 
specific way. Tensile lateral strain usually decreases r±2, while compression increases ryi- 
The bonding-relevant d xz and d yz orbitals of the top layer have the strongest interaction 
with the second layer atoms, so the relaxation of rvi significantly reduces the effect of lateral 
strain for these orbitals. This is why de xzyz /ds is less than ded/ds for all surfaces studied 
(Table P). 

The effect of r 12 relaxation on strain tunability is also strongly facet-dependent. On the 
more open (100) surface, relaxations of ri2 are larger, making de xzyz smaller for each (100) 
facet studied than for the corresponding (111) facet. This explains why the tunability of 
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E c hem (fix and rlx) and E^ ssoc are much lower on (100) surfaces than on (111) surfaces. 

Our demonstration that a single orbital-specific chemisorption model can be applied to 
different facets, strains, and metals, has implications for the modeling and design of more 
realistic catalyst surfaces. DFT studies have found that reactions on late TMs are more 
likely to proceed on defects such as steps and kinks. The model presented suggests that 
one should examine how different nearest-neighbor and inter-planar separations affect the 
orbital-specific electronic structure, and predict chemisorption properties accordingly. 

Incorporating the effects of strain and r 12 relaxation on the relevant PDOS centers greatly 
improved chemisorption modeling. We therefore suggest that if further couplings between 
PDOS centers and adsorbate structure can be parameterized, the resulting model could offer 
even greater accuracy and broader applicability. 

In conclusion, we use the energy levels of the substrate d-band projected onto the sub- 
strate atomic orbitals and their overlap with the CO bonding molecular orbitals, in a second- 
order perturbation theory-type model for chemisorption. The resulting model is able to ac- 
count for changes in -Ediem on different surface facets under different conditions of strain. We 
have also shown that trends in the dissociative chemisorption of CO at bridge site are gov- 
erned by the same orbital-specific factors. The results shown here should be generally valid 
for other molecular and atomic adsorption on higher index surfaces and for perturbations 
other than strain. 
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FIG. 1: (a) Variation in as a function of projection sphere cutoff radius r cut for Pt(lll). e x 2_ y 2, 
e xz , e z 2, e yz , and e xy PDOS centers are shown by circles, squares, diamonds, up-triangles, and 
down-triangles, respectively, (b) Same as (a) for Pt(100). 

The graphs demonstrate the importance of extrapolating r cu t — > a.u.: for Pt(lll) the asymptotic 
(i-band centers are more nearly equal than their large-r cut estimates; for Pt(100) they are more 
dissimilar. 
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FIG. 2: Plots of and E^? soc vs. and e xzyz for Pt (circle), Rh (square), and Pd (diamond) 

(111) surfaces (open) and (100) surfaces (filled). Data for five lateral strain states (0%, ±1%, and 
±2%) are shown. Linear regressions are shown for each metal, (a) -E^hcm vs. e<j, r cut =2 a.u. (b) 
£fssoc vs. e d , r cut =2 a.u. (c) E^ cm vs. e xzyz , r cut -0 a.u. (d) £^ soc vs. e d and e xzyz , r cut a.u. 
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FIG. 3: Correlation plots of modeled £" c hem and -E'chem.DFT- Pt (circle), Rh (square), and Pd 
(diamond) (111) surfaces (open) and (100) surfaces (filled). Data for five lateral strain states (0%, 
±1%, and ±2%) are shown. (a)^ m (Equation EJ) vs. E r c £ m DFT . (b)Eg s em (Equation [TJ vs. 

rprlx 

-^chem.DFT- 

The plots show that orbital-specific modeling (Equation ^) is required to achieve the same quality 
of correlation for both (111) and (100) facets. 
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